{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PHCv2.4.86 released 2022-05-30 works!\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pickle\n",
    "from scipy import io as sio\n",
    "import c_solver_2D as csolver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy import spatial\n",
    "from scipy import linalg as la"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using matplotlib backend: QtAgg\n"
     ]
    }
   ],
   "source": [
    "%matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "import BerryPy.TightBinding as tb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2D topology"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data/GBZ_data_6/x_major_1.pkl\", \"rb\") as fp:\n",
    "    all_segs = pickle.load(fp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "os.mkdir(\"data/GBZ_data_6_filtered\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "for j in range(50):\n",
    "    with open(\"data/GBZ_data_6/x_major_%d.pkl\"%(j), \"rb\") as fp:\n",
    "        all_segs = pickle.load(fp)\n",
    "    all_long_segs = []\n",
    "    for seg in all_segs:\n",
    "        if(seg.shape[0] > 1):\n",
    "            all_long_segs.append(seg)\n",
    "    with open(\"data/GBZ_data_6/x_major_%d_filtered.pkl\"%(j), \"wb\") as fp:\n",
    "        pickle.dump(all_long_segs, fp)\n",
    "    with open(\"data/GBZ_data_6/y_major_%d.pkl\"%(j), \"rb\") as fp:\n",
    "        all_segs = pickle.load(fp)\n",
    "    all_long_segs = []\n",
    "    for seg in all_segs:\n",
    "        if(seg.shape[0] > 1):\n",
    "            all_long_segs.append(seg)\n",
    "    with open(\"data/GBZ_data_6/y_major_%d_filtered.pkl\"%(j), \"wb\") as fp:\n",
    "        pickle.dump(all_long_segs, fp)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.gca()\n",
    "N_total_segs = 0\n",
    "upper_band = []\n",
    "lower_band = []\n",
    "for j in range(50):\n",
    "    with open(\"data/GBZ_data_6/x_major_%d_filtered.pkl\"%(j), \"rb\") as fp:\n",
    "        all_segs = pickle.load(fp)\n",
    "    for seg in all_segs:\n",
    "        N_total_segs += 1\n",
    "        ax.plot(seg[:,0].real, seg[:,0].imag)\n",
    "        if(seg[0,0].real > 0):\n",
    "            upper_band.append(seg)\n",
    "        else:\n",
    "            lower_band.append(seg)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "13404\n",
      "6676\n",
      "6728\n"
     ]
    }
   ],
   "source": [
    "print(N_total_segs)\n",
    "print(len(upper_band))\n",
    "print(len(lower_band))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "for seg in upper_band:\n",
    "    plt.plot(seg[:,0].real, seg[:,0].imag, 'r.')\n",
    "for seg in lower_band:\n",
    "    plt.plot(seg[:,0].real, seg[:,0].imag, 'b.')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "upper_band_arr = np.vstack(upper_band)\n",
    "lower_band_arr = np.vstack(lower_band)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 264,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f9b31694d60>]"
      ]
     },
     "execution_count": 264,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "OBC_data = sio.loadmat(\"data/OBC_spectrum_6_N150_N250.mat\")\n",
    "eigv = OBC_data['eigv'].flatten()\n",
    "plt.plot(eigv.real, eigv.imag, '.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 266,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(upper_band_arr[:,0].real, upper_band_arr[:,0].imag, 'b.')\n",
    "plt.plot(lower_band_arr[:,0].real, lower_band_arr[:,0].imag, 'r.')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "upper_band_k1 = np.log(upper_band_arr[:,1]).imag\n",
    "upper_band_k2 = np.log(upper_band_arr[:,2]).imag\n",
    "lower_band_k1 = np.log(lower_band_arr[:,1]).imag\n",
    "lower_band_k2 = np.log(lower_band_arr[:,2]).imag\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[2 3]\n",
      " [4 5]\n",
      " [6 7]\n",
      " [8 9]]\n"
     ]
    }
   ],
   "source": [
    "A = np.array([\n",
    "    [1,2],\n",
    "    [3,4],\n",
    "    [5,6],\n",
    "    [7,8]\n",
    "])\n",
    "print(A + np.array([1,1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 241,
   "metadata": {},
   "outputs": [],
   "source": [
    "def toric_Delaunay(points, X_period, Y_period):\n",
    "    N_points = points.shape[0]\n",
    "    # shifts\n",
    "    enlarged_points = np.vstack([\n",
    "        points,\n",
    "        points + np.array([X_period,0]),\n",
    "        points + np.array([0, Y_period]),\n",
    "        points + np.array([X_period,Y_period]),\n",
    "        points + np.array([-X_period,0]),\n",
    "        points + np.array([0, -Y_period]),\n",
    "        points + np.array([-X_period,Y_period]),\n",
    "        points + np.array([X_period, -Y_period]),\n",
    "        points + np.array([-X_period, -Y_period])\n",
    "    ])\n",
    "    enlarged_tri = spatial.Delaunay(enlarged_points)\n",
    "\n",
    "    # plt.triplot(enlarged_points[:,0], enlarged_points[:,1], enlarged_tri.simplices)\n",
    "    # plt.plot([-X_period/2, -X_period/2, X_period/2, X_period/2, -X_period/2], \n",
    "    #          [-Y_period/2, Y_period/2, Y_period/2, -Y_period/2, -Y_period/2], 'k-')\n",
    "\n",
    "    toric_simplices = []\n",
    "    # toric_simplices_extended = []\n",
    "    for curr_simplex in enlarged_tri.simplices:\n",
    "        curr_regions = set([j//N_points for j in curr_simplex])\n",
    "        cases = [\n",
    "            (0 in curr_regions) and (0 in curr_regions),\n",
    "            (0 in curr_regions) and (1 in curr_regions),\n",
    "            (0 in curr_regions) and (2 in curr_regions),\n",
    "            (0 in curr_regions) and (3 in curr_regions),\n",
    "            (1 in curr_regions) and (2 in curr_regions)\n",
    "        ]\n",
    "        if(any(cases) and curr_regions.issubset({0,1,2,3})):\n",
    "            curr_normalized_simplex = [(curr_index % N_points) for curr_index in curr_simplex]\n",
    "            toric_simplices.append(curr_normalized_simplex)\n",
    "            # toric_simplices_extended.append(list(curr_simplex))\n",
    "\n",
    "    #     if(any(np.asarray(curr_simplex) < N_points)):\n",
    "    #         if(not all(np.asarray(curr_simplex) < N_points)):\n",
    "    #             curr_normalized_simplex = [(curr_index % N_points) for curr_index in curr_simplex]\n",
    "    #             if(simplices_lookup.count(set(curr_normalized_simplex)) == 0):\n",
    "    #                 toric_simplices.append(curr_normalized_simplex)\n",
    "    #                 simplices_lookup.append(set(curr_normalized_simplex))\n",
    "    #                 toric_simplices_extended.append(list(curr_simplex))\n",
    "    #         else:\n",
    "    #             toric_simplices.append(list(curr_simplex))\n",
    "    #             simplices_lookup.append(set(curr_simplex))\n",
    "    #             toric_simplices_extended.append(list(curr_simplex))\n",
    "    # plt.triplot(enlarged_points[:,0], enlarged_points[:,1], toric_simplices_extended)\n",
    "    # plt.show()\n",
    "    return toric_simplices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f9b2f9975e0>]"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.plot(upper_band_k1, upper_band_k2, '.')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {},
   "outputs": [],
   "source": [
    "upper_band_tri = spatial.Delaunay(np.column_stack([upper_band_k1, upper_band_k2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 240,
   "metadata": {},
   "outputs": [],
   "source": [
    "upper_band_toric = toric_Delaunay(np.column_stack([upper_band_k1, upper_band_k2]), 2 * np.pi, 2 * np.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 244,
   "metadata": {},
   "outputs": [],
   "source": [
    "lower_band_toric = toric_Delaunay(np.column_stack([lower_band_k1, lower_band_k2]), 2 * np.pi, 2 * np.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "for curr_simplex in upper_band_tri.simplices:\n",
    "    curr_det = np.array(\n",
    "        [\n",
    "            [upper_band_k1[curr_simplex[1]] - upper_band_k1[curr_simplex[0]], upper_band_k1[curr_simplex[2]] - upper_band_k1[curr_simplex[0]]],\n",
    "            [upper_band_k2[curr_simplex[1]] - upper_band_k2[curr_simplex[0]], upper_band_k2[curr_simplex[2]] - upper_band_k2[curr_simplex[0]]]\n",
    "        ]\n",
    "    )\n",
    "    if(la.det(curr_det) < 0):\n",
    "        print(la.det(curr_det))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(projection='3d')\n",
    "# ax = fig.gca()\n",
    "# ax.plot_trisurf(upper_band_k1, upper_band_k2, upper_band_arr[:,0].real, triangles = upper_band_tri.simplices)\n",
    "ax.plot_trisurf(upper_band_k1, upper_band_k2, upper_band_arr[:,0].real, triangles = upper_band_toric)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 268,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 268,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(projection='3d')\n",
    "b = 5\n",
    "a = 1\n",
    "X = (b + a * np.cos(upper_band_k2)) * np.cos(upper_band_k1)\n",
    "Y = (b + a * np.cos(upper_band_k2)) * np.sin(upper_band_k1)\n",
    "Z = a * np.sin(upper_band_k2)\n",
    "# ax.plot(upper_band_k1, upper_band_k2, upper_band_arr[:,0].real, '.')\n",
    "# ax.plot_trisurf(X, Y, Z, triangles = upper_band_tri.simplices)\n",
    "ax.plot_trisurf(X, Y, Z, triangles = upper_band_toric, color=(0,0,0,0), edgecolor=(0,0,0,1))\n",
    "ax.set_aspect(\"equal\")\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "ax.set_zticks([])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data/GBZ_data_7/x_major_30.pkl\", \"rb\") as fp:\n",
    "    all_segs = pickle.load(fp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(projection='3d')\n",
    "for seg in all_segs:\n",
    "    ax.plot(seg[:,1].real, seg[:,1].imag, seg[:,0].real)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1204\n"
     ]
    }
   ],
   "source": [
    "print(len(all_segs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "one_point_segs = []\n",
    "for seg in all_segs:\n",
    "    if(seg.shape[0] == 1):\n",
    "        one_point_segs.append(seg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "736\n"
     ]
    }
   ],
   "source": [
    "print(len(one_point_segs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1204\n"
     ]
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(projection='3d')\n",
    "print(len(all_segs))\n",
    "for seg in all_segs:\n",
    "    ax.plot(np.log(seg[:,1]).imag,np.log(seg[:,2]).imag , seg[:,0].real)\n",
    "for seg in one_point_segs:\n",
    "    ax.plot(np.log(seg[:,1]).imag,np.log(seg[:,2]).imag , seg[:,0].real, '.')\n",
    "    # ax.plot(seg[:,1].real, seg[:,1].imag, seg[:,0].real)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "loops, fragments = csolver.segments_to_loops(all_segs, 0.5, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "18\n",
      "240\n"
     ]
    }
   ],
   "source": [
    "print(len(loops))\n",
    "print(len(fragments))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "199\n",
      "400\n",
      "400\n",
      "203\n",
      "150\n",
      "307\n",
      "53\n",
      "174\n",
      "42\n",
      "24\n",
      "70\n",
      "18\n",
      "120\n",
      "36\n",
      "18\n",
      "136\n",
      "26\n",
      "24\n"
     ]
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(projection='3d')\n",
    "for curr_loop in loops:\n",
    "    # ax.plot(curr_loop[:,1].real, curr_loop[:,1].imag, curr_loop[:,0].real)\n",
    "    ax.plot(np.log(curr_loop[:,1]).imag,np.log(curr_loop[:,2]).imag , curr_loop[:,0].real)\n",
    "    print(curr_loop.shape[0])\n",
    "    # plt.pause(1)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n",
      "2\n",
      "2\n",
      "2\n"
     ]
    }
   ],
   "source": [
    "for curr_loop in loops:\n",
    "    print(csolver.winding_number(curr_loop[:,1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_eigenstates(model:tb.TightBindingModel, GBZ_point):\n",
    "    H = model.get_bulk_Hamiltonian_dense(np.log(GBZ_point[1:])/(2j*np.pi))\n",
    "    eigv_eq = np.eye(H.shape[0]) * GBZ_point[0] - H\n",
    "    eigvec = la.null_space(eigv_eq, rcond=1e-6)\n",
    "    # print(la.det(eigv_eq))\n",
    "    # print(eigvec)\n",
    "    return eigvec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [],
   "source": [
    "def winding_phase_fun(model:tb.TightBindingModel, vec_list:list):\n",
    "    winding_phase_multiplier = 1\n",
    "    # vec_list = [get_eigenstates(model, point_list[j,:])[:,0] for j in loop_indices]\n",
    "    for j in range(len(vec_list) - 1):\n",
    "        winding_phase_multiplier *= (vec_list[j].conjugate()).dot(vec_list[j+1])\n",
    "    winding_phase_multiplier *= (vec_list[-1].conjugate()).dot(vec_list[0])\n",
    "    return -np.log(winding_phase_multiplier).imag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Chern_number(model:tb.TightBindingModel, point_list:np.ndarray, simplices:list):\n",
    "    Berry_flux = 0\n",
    "    all_vecs = [get_eigenstates(model, point_list[j,:])[:,0] for j in range(point_list.shape[0])]\n",
    "    for curr_simplex in simplices:\n",
    "        curr_vec_list = [all_vecs[j] for j in curr_simplex]\n",
    "        Berry_flux += winding_phase_fun(model, curr_vec_list)\n",
    "    return Berry_flux / (2*np.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 195,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.08631383151736598\n"
     ]
    }
   ],
   "source": [
    "print(Chern_number(model, upper_band_arr, upper_band_tri.simplices))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 243,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.9999999999999962\n"
     ]
    }
   ],
   "source": [
    "print(Chern_number(model, upper_band_arr, upper_band_toric))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 245,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9999999999999972\n"
     ]
    }
   ],
   "source": [
    "print(Chern_number(model, lower_band_arr, lower_band_toric))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [],
   "source": [
    "def area_calculation(k1_list, k2_list, toric_simplices:list):\n",
    "    all_area = 0\n",
    "    for curr_simplex in toric_simplices:\n",
    "        p1x = k1_list[curr_simplex[1]] - k1_list[curr_simplex[0]]\n",
    "        p_list = [p1x - 2 * np.pi, p1x, p1x + 2 * np.pi]\n",
    "        min_p = np.argmin(np.abs(p_list))\n",
    "        p1x = p_list[min_p]\n",
    "\n",
    "        p1y = k2_list[curr_simplex[1]] - k2_list[curr_simplex[0]]\n",
    "        p_list = [p1y - 2 * np.pi, p1y, p1y + 2 * np.pi]\n",
    "        min_p = np.argmin(np.abs(p_list))\n",
    "        p1y = p_list[min_p]\n",
    "\n",
    "        p2x = k1_list[curr_simplex[2]] - k1_list[curr_simplex[0]]\n",
    "        p_list = [p2x - 2 * np.pi, p2x, p2x + 2 * np.pi]\n",
    "        min_p = np.argmin(np.abs(p_list))\n",
    "        p2x = p_list[min_p]\n",
    "\n",
    "        p2y = k2_list[curr_simplex[2]] - k2_list[curr_simplex[0]]\n",
    "        p_list = [p2y - 2 * np.pi, p2y, p2y + 2 * np.pi]\n",
    "        min_p = np.argmin(np.abs(p_list))\n",
    "        p2y = p_list[min_p]\n",
    "\n",
    "        all_area += la.det([\n",
    "            [p1x, p2x],\n",
    "            [p1y, p2y]\n",
    "        ])/2\n",
    "    return all_area"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 242,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "39.4784176043435\n"
     ]
    }
   ],
   "source": [
    "print(area_calculation(upper_band_k1, upper_band_k2, upper_band_toric))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f9b4c7733a0>]"
      ]
     },
     "execution_count": 119,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# plt.plot(upper_band_k1/np.pi, upper_band_k2/np.pi, '.')\n",
    "plt.triplot(np.column_stack([upper_band_k1/np.pi, upper_band_k2/np.pi]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "39.47841760435743\n"
     ]
    }
   ],
   "source": [
    "print((2*np.pi) **2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data/GBZ_model_6.pkl\", \"rb\") as fp:\n",
    "\n",
    "    model = pickle.load(fp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4.440172665480811e-17+6.5527436948497284e-18j)\n"
     ]
    }
   ],
   "source": [
    "eigvec = get_eigenstates(model, upper_band_arr[0,:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9999999999999998\n"
     ]
    }
   ],
   "source": [
    "print(la.norm(eigvec))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 269,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '$\\\\mathrm{Im}\\\\log(\\\\beta_2)$')"
      ]
     },
     "execution_count": 269,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N_points = 300\n",
    "k1_list = np.random.rand(N_points) * (2 * np.pi) - np.pi\n",
    "k2_list = np.random.rand(N_points) * (2 * np.pi) - np.pi\n",
    "plt.plot(k1_list, k2_list, '.')\n",
    "plt.xlabel(\"$\\\\mathrm{Im}\\\\log(\\\\beta_1)$\")\n",
    "plt.ylabel(\"$\\\\mathrm{Im}\\\\log(\\\\beta_2)$\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 270,
   "metadata": {},
   "outputs": [],
   "source": [
    "toric_test = toric_Delaunay(np.column_stack([k1_list, k2_list]), 2 * np.pi, 2 * np.pi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 275,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 275,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(projection='3d')\n",
    "b = 5\n",
    "a = 1\n",
    "X = (b + a * np.cos(k2_list)) * np.cos(k1_list)\n",
    "Y = (b + a * np.cos(k2_list)) * np.sin(k1_list)\n",
    "Z = a * np.sin(k2_list)\n",
    "# ax.plot(k1_list, k2_list, upper_band_arr[:,0].real, '.')\n",
    "# ax.plot_trisurf(X, Y, Z, triangles = upper_band_tri.simplices)\n",
    "ax.plot_trisurf(X, Y, Z, triangles = toric_test, color=(0,0,0,0), edgecolor=(0,0,0.5,1))\n",
    "ax.set_aspect('equal')\n",
    "ax.set_xticks([])\n",
    "ax.set_yticks([])\n",
    "ax.set_zticks([])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f9b510080a0>,\n",
       " <matplotlib.lines.Line2D at 0x7f9b1ef47a90>]"
      ]
     },
     "execution_count": 184,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plt.triplot(k1_list, k2_list, toric_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 188,
   "metadata": {},
   "outputs": [],
   "source": [
    "plane_test = spatial.Delaunay(np.column_stack([k1_list, k2_list]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 239,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "39.47841760435746\n"
     ]
    }
   ],
   "source": [
    "print(area_calculation(k1_list, k2_list, toric_test))\n",
    "# print(area_calculation(k1_list, k2_list, plane_test.simplices))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 216,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "39.47841760435743"
      ]
     },
     "execution_count": 216,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(2*np.pi) ** 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Try CGAL"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 197,
   "metadata": {},
   "outputs": [
    {
     "ename": "ModuleNotFoundError",
     "evalue": "No module named 'cgal'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mModuleNotFoundError\u001b[0m                       Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[197], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[39mimport\u001b[39;00m \u001b[39mcgal\u001b[39;00m\u001b[39m.\u001b[39;00m\u001b[39mperiodic_2d_triangulation\u001b[39;00m \u001b[39mas\u001b[39;00m \u001b[39mcgal_pt\u001b[39;00m\n",
      "\u001b[0;31mModuleNotFoundError\u001b[0m: No module named 'cgal'"
     ]
    }
   ],
   "source": [
    "import cgal.periodic_2d_triangulation as cgal_pt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Minimal winding"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 248,
   "metadata": {},
   "outputs": [],
   "source": [
    "with open(\"data/GBZ_data_8/y_major_0.pkl\", \"rb\") as fp:\n",
    "    all_segments = pickle.load(fp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 251,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = sio.loadmat(\"data/GBZ_char_poly_8.mat\")\n",
    "coeffs = data['coeffs'].flatten()\n",
    "degrees = data['orders'].flatten()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 253,
   "metadata": {},
   "outputs": [],
   "source": [
    "char_poly = csolver.CLaurant(3)\n",
    "char_poly.set_Laurant_by_terms(csolver.CScalarVec(coeffs), csolver.CIntVec(degrees))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 254,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+(-1+0*i)*(y**2)+(-1+0*i)*x*y+(-1+0*i)*(x**2)+(-1+0*i)*(x**2)*(y**2)+(-1+0*i)*(x**3)+(-1+0*i)*(x**3)*y+(1+0*i)*E*(x**2)*y;\n"
     ]
    }
   ],
   "source": [
    "print(char_poly.num.to_str(csolver.CStrVec([\"E\", \"x\", \"y\"])))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 255,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "aGBZ_2D_solver_IndexVec[0, 2, 1]\n"
     ]
    }
   ],
   "source": [
    "print(char_poly.denom_orders)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 256,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(projection='3d')\n",
    "for seg in all_segments:\n",
    "    ax.plot(seg[:,2].real, seg[:,2].imag, seg[:,0].real)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
